2′-Fucosyllactose helps butyrate producers outgrow competitors in infant gut microbiota simulations

Summary A reduced capacity for butyrate production by the early infant gut microbiota is associated with negative health effects, such as inflammation and the development of allergies. Here, we develop new hypotheses on the effect of the prebiotic galacto-oligosaccharides (GOS) or 2′-fucosyllactose (2′-FL) on butyrate production by the infant gut microbiota using a multiscale, spatiotemporal mathematical model of the infant gut. The model simulates a community of cross-feeding gut bacteria in metabolic detail. It represents the community as a grid of bacterial populations that exchange metabolites, using 20 different subspecies-specific metabolic networks taken from the AGORA database. The simulations predict that both GOS and 2′-FL promote the growth of Bifidobacterium, whereas butyrate producing bacteria are only consistently abundant in the presence of propane-1,2-diol, a product of 2′-FL metabolism. In absence of prebiotics or in presence of only GOS, however, Bacteroides vulgatus and Cutibacterium acnes outcompete butyrate producers by consuming intermediate metabolites.


Highlights
A computational model simulates bacterial ecology in the infant gut Model predicts that the prebiotics GOS and 2 0 -FL stimulate Bifidobacterium

Bacteroides vulgatus may outcompete butyrate producers
The prebiotic 2 0 -FL may help butyrate producers outcompete other species

INTRODUCTION
Infants develop a complex microbiota shortly after birth, which is important for healthy growth and development. 1Here, we focus on butyrate, a short-chain fatty acid (SCFA) that is produced in significant amounts by the gut bacteria 2 and is absorbed by the gut colonocytes.Production of butyrate by the microbiota has been suggested to improve the health of infants in a number of ways.Firstly, butyrate in the gut is a key energy source for the gut epithelium, making it important for maintaining the gut barrier function. 3A breakdown of the gut barrier function due to a lack of butyrate is associated with diseases such as inflammatory bowel disease and rectal cancer. 3,46][7] Infant butyrate producing bacteria provide protection against food allergies when transplanted into a mouse model, 8 suggesting causality.Butyrate production is also associated with a reduced risk of colic in infants. 9Butyrate also modulates the immune system throughout the body, inhibiting inflammation and carcinogenesis. 10These data suggest it may be desirable to stimulate butyrate production in the infant gut.Using mechanistic computational modeling, here we investigate how stimulation of butyrate producing bacteria may be achieved in the early infant gut microbiota through supplementation with prebiotics.
Microbiota composition and metabolism are influenced by endogenous factors, e.g., gut maturity and inflammation, and exogenous factors, e.g., nutrition, probiotics, and antibiotics.Here, we focus on nutrition, which is the primary exogenous factor.3][14] It has been hypothesized that some of the health effects associated with prebiotics may be linked to indirect stimulation of butyrate producing bacteria. 7,15][14] Butyrate producing bacteria such as Anaerobutyricum hallii (formerly Eubacterium hallii 16 ) cannot directly consume GOS or 2 0 -FL, but they can consume metabolites of GOS or 2 0 -FL digestion. 17The primary consumers of GOS and 2 0 -FL in the infant gut are Bifidobacterium spp. 18,19etabolites produced by Bifidobacterium spp., in turn, become important food sources for butyrate producing bacteria.For example, in vitro it has been found that the butyrate producing bacterium A. hallii can feed on lactate and propane-1,2-diol (1,2-PD), which are metabolites of Bifidobacterium spp. 17A. hallii can also coexist with Bifidobacterium longum ssp.infantis in vitro on a substrate of glucose or 2 0 -FL. 17espite these in vitro findings that demonstrate potential coexistence of Bifidobacterium spp.and butyrate producing bacteria, in vivo, i.e., in the infant gut microbiota, butyrate producing bacteria often only have a low abundance and butyrate is found in the feces of only 35% of infants. 20It is unclear why butyrate producing bacteria and butyrate are not commonly abundant in vivo, given that in vitro cross-feeding on lactate occurs readily, 17 and that lactate-producing Bifidobacterium species are abundant in the gut of most infants. 21,22Using computational modeling, we explore the conditions that may stimulate butyrate producing bacteria in vivo in the infant gut.To this end, we will compare simulations of simple microbial communities, such as those studied in vitro, with simulations of more complex communities that may more closely resemble the in vivo situation.
Briefly, the computational model suggests that in simple microbial communities, populations of butyrate producing bacteria can crossfeed on Bifidobacterium metabolites.However, in more complex communities the intermediary metabolites are consumed by competitors instead of butyrate producing bacteria.In the presence of 2 0 -FL, populations of butyrate producing bacteria are nevertheless supported.The mechanism suggested by our simulations is that Bifidobacterium produces 1,2-PD from 2 0 -FL, which specifically feeds butyrate producing species, allowing these to outgrow competing cross-feeders.We provide predictions for interactions in in vivo and in vitro systems and suggestions for in vitro verification of these predictions.

Model outline
To develop new hypotheses on how oligosaccharides can stimulate the production of butyrate, we further develop a multiscale metabolic model (Figures 1A and 1B) of the carbon metabolism of the infant gut microbiota. 23The computational model is based on our earlier models of the adult and infant microbiota. 23,24In comparison with these previous models, the present model simulates a larger number of small bacterial populations, using a larger, more diverse, and further curated set of metabolic models of gut bacteria from the AGORA database. 25In particular, we have included the butyrate producers A. hallii, Roseburia inulinivorans, and Clostridium butyricum and the digestion of the prebiotic oligosaccharides GOS and 2 0 -FL by Bifidobacterium longum ssp.infantis.The complete community model integrates these predictions of metabolism over space and time to create a multiscale model that covers the development and variation of the infant gut microbiota over the first three weeks of life.Other multiscale metabolic modeling techniques have been used previously to model the adult human microbiota in frameworks such as SteadyCom and Comets. 26,27The model presented here distinguishes itself from these frameworks by its focus on the infant gut microbiota, by including factors such as prebiotics and the initial presence of oxygen at birth.
Briefly, the spatial model simulates the ecology of an intestinal microbial ecosystem and features genome-scale metabolic models (GEMs) of intestinal bacteria, spatial structuring, exchange of extracellular metabolites, and population dynamics.The system is simulated on a regular square lattice of 22538 boxes of 232 mm, representing a typical infant colon of 4531:6 cm.Each box contains a simulated metapopulation of one of a set of up to 20 of the most common bacterial species present in the infant gut 21 (Table 1), and concentrations of simulated nutrients and metabolites such as extracellular oligosaccharides and short-chain fatty acids.Based on the concentrations of metabolites, the systems predict the growth rate for each metapopulation as well as the uptake and excretion rates of metabolites using a GEM taken from AGORA, 28 a database of metabolic networks of intestinal bacteria.The system is initialized by distributing, on average, 540 populations over the system at random.Oxygen is introduced during initialization, and water is always available.
After initialization, the model is simulated in timesteps representing 3 min of real time.Each timestep of the simulation proceeds as follows.Every 3 h (i.e., 60 timesteps), a mixture of simulated lactose and/or oligosaccharides is added to the leftmost six columns of lattice sites.Then, in each step, the model predicts the metabolism of each local population using flux balance analysis (FBA) based on the metabolites present in the local lattice site, the GEM of the species, and the enzymatic constraint.The enzymatic constraint limits the total amount of metabolism that can be performed by each local population per timestep by limiting the maximum summed flux for each FBA solution.The enzymatic constraint scales linearly with the local population size.This approach allows us to model metabolic switches and trade-offs. 23,29The FBA solution includes a set of influx rates and efflux rates of metabolites that are used to update the environmental metabolite concentrations.The local populations are assumed to grow at a rate linearly proportional to the rate of ATP production, 30 which is predicted by FBA by optimizing for ATP production rates.Populations may create a new population in a neighboring lattice site if the local population is 200 times the initial size (Figure 1A-1).Populations of more than 400 times the local size, which can only form when density if so high new populations cannot be created, stop metabolism to represent quiescence.Populations spread at random into adjacent lattice sites (Figure 1A-2); metabolites diffuse and advect toward the back of the tube (Figure 1A-3&4).To mimic excretion, metabolites and populations are deleted from the most distal column each timestep.To represent bacterial colonization, new populations of randomly selected species are introduced into empty lattice sites at a small probability.All parameters are given in Table 2. Details of the model are given in the STAR Methods.

Model with simplified consortium of species predicts coexistence of butyrate producing bacteria and Bifidobacterium
We first simulated the model using a simplified consortium of species, the two Bifidobacterium longum subspecies (Table 1) and three butyrate producing species: Anaerobutyricum hallii, Clostridium butyricum, and Roseburia inulinivorans.We performed 30 simulations for each of four conditions, in which the following sugars were introduced every three simulated hours: (1) 211 mmol lactose and no prebiotics, (2) 422 mmol lactose and no prebiotics, (3) 211 mmol lactose plus 211 mmol GOS, and (4) 211 mmol lactose plus 211 mmol 2 0 -FL.We estimated 211 mmol lactose to be a realistic amount of lactose to reach the infant colon, given infant intake and small intestinal uptake. 26,31As there is little absorption by the small intestine of prebiotics, 32 the amount of prebiotics in the nutrition consumed by the infant would be much smaller than the amount of lactose.We also include the 422 mmol lactose condition to control for the possibility that effects in the conditions with prebiotics are due the larger amount of sugar present conditions, instead of their type.The condition with 422 mmol lactose does not correspond to an in vivo condition.We analyzed the abundance of each species at the end of 10,080 timesteps, representing 21 simulated days.In each of the four conditions, Bifidobacterium bacteria (Figure 1C) and butyrate producing bacteria coexisted (Figure 1D), and, paradoxically, butyrate producing bacteria were reduced in presence of prebiotics.
In the presence of competitors, model predicts coexistence of butyrate producing bacteria and Bifidobacterium in the presence of 2'-FL but not in the presence of GOS We next examined the behavior of the system in the presence of a more complex consortium, consisting of all 20 species and subspecies listed in Table 1, simulating the same four conditions.In the absence of prebiotics, regardless of the quantity of lactose, the model predicted that Bifidobacterium, Bacteroides, and Escherichia became the most abundant genera after three weeks (Figure 2A, Video S1), consistent with in vivo observation. 21,22We also observed some abundance of Bacilli in accordance with in vivo observations. 21,22,33The higher quantity of (C and D) Abundance of (C) Bifidobacterium spp., (D) butyrate producing bacteria, at the end of 21 days for 30 sets of simulations with no prebiotics, no prebiotics and additional lactose, with GOS, or with 2'-FL at the end of 21 days n = 30 for each condition, each simulation is represented by one dot.See also Table S1.lactose resulted in a higher average abundance for all major groups.In the absence of prebiotics, butyrate producing bacteria achieved a combined abundance over 1$10 10 in only 4 of the 30 simulations with 211 mmol of lactose per 3 h, and 6 of the 30 with 422 mmol of lactose (Figure 2B).In the remaining simulations, the butyrate producing bacteria remained almost absent, staying below 1$10 10 bacteria.In the simulations with GOS, Bifidobacterium was more abundant than in the condition without prebiotics (p < 0.001, Figure 2A) whereas the butyrate producing bacteria were not affected (p = 0.18) (Figure 2B).With GOS, butyrate producing bacteria also had a combined abundance of over 1$10 10 bacteria at the end of 13 of the 30 simulations (Figure 2B).Interestingly, in the condition with 2 0 -FL the abundance of butyrate producing bacteria was over 1$10 10 bacteria at the end of 19 of 30 simulations (Figure 2B), and the butyrate producing species were more abundant (Figure 2A, Video S2) than in the other conditions.Thus 2 0 -FL but not GOS stimulated butyrate producing bacteria in the complex community.To test for any concentration dependence or crosstalk between 2 0 -FL and GOS, we next performed sets of 30 simulations in the presence of 211 mmol lactose and levels of 2 0 -FL and GOS varying between 21.1 mmol and 211 mmol per 3 h and combinations thereof (Figure S1).The amount of 2 0 -FL (p = 0.017, Kruskal-Wallis rank-sum test) but not that of GOS (p = 0.658, Kruskal-Wallis rank-sum test) affected the abundance of butyrate producing bacteria, further supporting the prediction that 2 0 -FL but not GOS stimulates butyrate producing bacteria in the complex community.
In order to investigate why 2 0 -FL led to a more consistent abundance of butyrate producing bacteria, we analyzed the metabolic interactions between bacterial species.We visualized the network of metabolic fluxes between the bacteria using arrows between species and metabolite pools in Figures 2C-2E.The resulting diagrams show both primary consumption, i.e., uptake of nutrients such as lactose, GOS, and 2 0 -FL, and cross-feeding, i.e., uptake of metabolites produced by other species.Sample visualizations of the condition without prebiotics (211 mmol lactose) (Figure 2C, Video S3) and the condition with GOS (Figure 2D) revealed co-occurrence of species and cross-feeding but no butyrate production.In these simulations, the cross-feeding metabolite lactate, which is a known substrate for butyrate producing bacteria, 17 was consumed by Bacteroides vulgatus and Cutibacterium acnes, respectively.Butyrate formation only occurred in the sample simulation with 2 0 -FL (Figure 2E).Only in the presence of 2 0 -FL and not in the other conditions, was a flux of 1,2-PD directed toward the butyrate producing species (Figure 2E; Video S4).We, therefore, hypothesized that butyrate producing species may be more abundant in the model simulations with 2 0 -FL because 2 0 -FL digestion by Bifidobacterium produces 1,2-PD as a cross-feeding substrate.1,2-PD is a known Bifidobacterium metabolite from 2 0 -FL in vitro. 17To test this hypothesis, we performed new sets of simulations with 2 0 -FL in which we blocked the uptake by butyrate producing bacteria of either lactose, lactate, or 1,2-PD, i.e., the uptake of metabolites most consumed by butyrate producing bacteria was disabled.Indeed, blocking the uptake of any of these metabolites led to a reduction of butyrate producing bacteria (Figure 2F).Thus a flux of not only lactose and lactate, but also 1,2-PD that is only produced in presence of 2 0 -FL, was required for sustaining butyrate producing bacteria in our simulations.
We next turned to the model with the simplified consortium of species, the two Bifidobacterium subspecies and three butyrate producing species, to test if uptake of lactose, lactate, and 1,2-PD was also required for butyrate producing bacteria to become abundant with this consortium.After blocking the uptake of lactose, lactate, or 1,2-PD by butyrate producing bacteria, the abundance of butyrate producing bacteria was reduced at the end of the simulations compared to the control (Figure 3A).Surprisingly, however, and in contrast to the complete system (Figure 2F), butyrate producing populations retained an abundance of over 1$10 10 bacteria in respectively 27 and 30 of 30 simulations when lactose or 1,2-PD uptake was disabled.Thus neither lactose nor 1,2-PD were essential for butyrate producing bacteria.Altogether, 1,2-PD, and thus 2 0 -FL, was required for butyrate producing bacteria in the complex system but not in the simplified system.Thus, these model simulations suggest that supplementation with 2 0 -FL introduces a flux of 1,2-PD from Bifidobacterium spp. to butyrate producing bacteria that prevents competitive exclusion of butyrate producers by competitors such as B. vulgatus (Figure 2C) or C. acnes (Figure 2D).Bacteroides vulgatus and C. acnes are effective competitors on different substrates In the 2 0 -FL condition, butyrate producing bacteria fed on lactate and 1,2-PD (Figure 2E).In the conditions without 2 0 -FL, no 1,2-PD was produced and lactate was consumed by B. vulgatus or C. acnes (Figures 2C and 2D).This suggests that, in the absence of 1,2-PD, B. vulgatus and C. acnes outcompete the butyrate producing bacteria for lactate.To investigate whether these species could indeed be responsible for outcompeting butyrate producing bacteria, we again turned to the model with the simplified consortium and added the potential competitors B. vulgatus and C. acnes to the consortium one by one.First, we studied the simplified consortium in absence of prebiotics in the conditions with 211 mmol and 422 mmol lactose per 3 h.The abundance of butyrate producing bacteria was reduced in presence of B. vulgatus but not in presence of C. acnes (Figure 3B, 422 mmol visualized in Figure S3).After blocking lactose or lactate uptake by B. vulgatus in the condition with 211 mmol lactose, the abundance of butyrate producing bacteria was restored (Figure 3B), indicating that B. vulgatus required both lactose and lactate to effectively outcompete the butyrate producing bacteria.
In the conditions with GOS, the situation was reversed: C. acnes but not B. vulgatus outcompeted butyrate producing bacteria (Figure 3C).After blocking uptake of lactate by C. acnes the abundance of butyrate producing bacteria was restored (Figure 3C).C. acnes does not use lactose in the model.Taken together, these simulations suggest that lactate is required for competitive exclusion of butyrate producing bacteria by C. acnes.
In the condition with 2 0 -FL, B. vulgatus did not outcompete butyrate producing bacteria (Figure 3D).C. acnes (p = 0.001) moderately suppressed butyrate producing bacteria, with 29 of 30 simulations still predicting an abundance of butyrate producing bacteria of over 1$ 10 10 bacteria.This agrees with the simulations using the full consortium (Figure 2B), which also displayed a robust abundance of butyrate producing bacteria in the 2 0 -FL condition.
Butyrate producing bacteria can use a mixture of lactate and 1,2-PD as substrates in the 2 0 -FL condition to grow faster than their competitors To analyze how butyrate producing bacteria can outcompete other species only in the presence of 2 0 -FL but not in the presence of GOS or without prebiotics, we next examined the growth rates per timestep on unlimited quantities of the three key substrates for butyrate producing bacteria indicated previously: lactose, lactate, and 1,2-PD.With unlimited availability of lactose, the growth of the three butyrate producing species was reduced relative to the growth of most other species (Figure 4A).With unlimited lactate, growth for butyrate producing species was superior to the other species but not to C. acnes (Figure 4B).In presence of unlimited 1,2-PD and acetate the butyrate producing species A. hallii and Roseburia inulinivorans grew faster than the other species (Figure 4C).On a mixture of limited lactate and 1,2-PD, with acetate available, two of the three butyrate producing species also grew faster compared to all other species (Figure 4D).Thus the unique ability of butyrate producing bacteria to grow on 1,2-PD and acetate in the model allowed them to outcompete other lactate-consuming species in environments with 1,2-PD, such as those where Bifidobacterium consumes 2 0 -FL.However, they would be unable to outcompete the same species in conditions without 1,2-PD.

Sensitivity analysis
Finally, to test the generality of our observations, we performed a sensitivity analysis on the system.The enzymatic constraint (Figures S2A and  S2B), the death rate and growth rate (through the ATP required per population unit) (Figures S2C and S2D), the placement of new populations of random species in empty lattice sites (colonization) (Figures S2E and S2F), the diffusion of metabolites and populations (Figures S2G and  S2H), the amount of initial oxygen (Figures S2I and S2J), the presence of quiescence for large populations (Figure S2K), well-mixed conditions (Figure S2L), and the size of the lattices (Figures S2M and S2N) were varied.The exact implementation of colonization, well-mixed conditions, quiescence, and lattice size are described in the methods.We used three conditions for most changed parameters: 211 mmol lactose, 211 mmol lactose plus 211 mmol GOS, and 211 mmol lactose plus 211 mmol 2 0 -FL per 3 h.We only used the latter two for disabling quiescence, as no populations entered quiescence during our initial runs with 211 mmol lactose.We found minor sensitivity for most parameter changes (Figure S2).We found the most notable effects when we disabled colonization, or initial oxygen, when we used well-mixed conditions, and when we decreased the lattice size.When we disabled colonization, the abundance of butyrate producing bacteria was lower in all three conditions (p < 0.001 for all, Figure S2E).The absence of initial oxygen increased the abundance of butyrate producing bacteria in the condition without prebiotics and with 2 0 -FL (p = 0.002,p = 0.035, Figure S2I).When we used well-mixed conditions, which disabled spatial separation, Bifidobacterium spp.were much less abundant without prebiotics, and B. vulgatus was much more abundant (p < 0.001, Figure S2L).There were also larger butyrate producer populations with GOS (p = 0.012, Figure S2L) and with 2 0 -FL (p < 0.001, Figure S2L).This indicates that the butyrate producers were more competitive in this well-mixed environment in the model.With a smaller lattice size, each lattice site represented 434 mm of space, instead of 232 mm.In this condition, Bifidobacterium was much less abundant in all conditions (p < 0.001, Figure S2M), which represents a mismatch with in vivo data. 21,22Conversely, species in the ''others'' category were much more abundant (p < 0.001, Figure S2M), as were butyrate producers without (p = 0.01) and with prebiotics (p < 0.001).A larger lattice, where each lattice site represented 131 mm of space, led to very similar results as with our default lattice size (Figure S2N) but was much more computationally intensive.Taken together, these results indicate that sustained colonization, the presence of initial oxygen, non-well-mixed conditions, and a minimum lattice size are particularly important in the simulated system.

DISCUSSION
This paper describes a computational study of the effects of the prebiotics GOS and 2 0 -FL on butyrate producing bacteria in the infant gut microbiota.We have used the model to generate novel hypotheses to explain the-sometimes counter-intuitive-mechanisms at the biochemical and population level that underlie the effects of prebiotics.The model predicts that butyrate producing bacteria can coexist with Bifidobacterium in the infant gut with or without GOS or 2 0 -FL as long as no other bacterial species are present.As soon as other bacterial species are introduced into the model, we found that they can act as competitors, thus reducing the abundance of butyrate producing bacteria.Specifically, the model predicts that B. vulgatus outcompetes butyrate producing bacteria in absence of prebiotics.The predicted mechanism is that B. vulgatus consumes lactose and lactate, important food sources of the butyrate producing species.In presence of GOS, the model predicts that C. acnes becomes the key competitor of the butyrate producing bacteria due to its lactate consumption.In presence of 2 0 -FL, however, the butyrate producing species are no longer outcompeted.The mechanism as predicted by the model is as follows.The breakdown of 2 0 -FL by Bifidobacterium produces 1,2-PD.1,2-PD becomes an additional food source for the butyrate producing bacteria, helping them to outgrow competitors.Thus, our modeling results predict that only 2 0 -FL, but not GOS, supports populations of butyrate producing bacteria in their competition against species such as B. vulgatus and C. acnes.
The following in vitro and in vivo observations agree with these model predictions.Firstly, the model predicts coexistence and crossfeeding between Bifidobacterium and butyrate producing species on 2 0 -FL.In agreement with the model predictions, coexistence of and cross-feeding between Bifidobacterium and butyrate producing bacteria occurs in vitro within simplified, synthetic communities on glucose, fucose, and 2 0 -FL in the absence of competitors. 17Secondly, the model predicts that in presence of the competitors such as B. vulgatus and C. acnes, B. vulgatus will become abundant in the absence of prebiotics and outcompete butyrate producing species.In agreement with this model prediction, B. vulgatus is often abundant in the in vivo infant gut microbiota, 21 and it can consume lactose in vitro. 34No information is available on lactate consumption of B. vulgatus, but the related Bacteroides fragilis is able to consume lactate in vitro. 35Thirdly, the model predicts that C. acnes outcompetes butyrate producing bacteria in presence of GOS by consuming lactate.In agreement with this prediction, C. acnes is found in 22% of infants 21 and Cutibacterium avidum, closely related to C. acnes, 36 reduces the abundance of the butyrate producer A. hallii in an in vitro lactate-fed microbiota from infant fecal samples. 37Both C. acnes and C. avidum consume lactate in vitro. 38Finally, the model predicts that butyrate producing bacteria become competitive through cross-feeding on 1,2-PD, which is produced by Bifidobacterium longum from 2 0 -FL.In agreement with this prediction, the butyrate producer A. hallii cross-feeds on 1,2-PD in an in vitro synthetic community of A. hallii and B. longum. 17Also in line with this prediction, 2 0 -FL supplementation increased the abundance of butyrate producing bacteria in in vitro fecal communities based on infant fecal samples, which likely include key competitors of butyrate producing species. 15An in vitro colonic fermentation model inoculated with infant feces has previously been used to study the effect of introducing specific competitors to lactate-consuming infant gut microbiota. 37This approach could also be used to test if B. vulgatus and C. acnes are viable competitors in the infant gut and if the presence of 1,2-PD allows butyrate producing species to outcompete other bacteria.
More broadly, the model simulations without prebiotics predict that Escherichia, Bacteroides, and Bifidobacterium become the three most abundant genera, which agrees with the most abundant genera found in the infant gut microbiota around the age of three weeks. 21,22The relative abundances the model predicts for butyrate producing species range from 1.4% without prebiotics to 4.8% with 2 0 -FL, both of which are within the broad range of values reported for the butyrate producing community. 20However, for two less abundant groups, Bacilli and Veillonella, the model predictions disagree with in vivo data.Firstly, an initially dominant Bacilli phase is sometimes seen in vivo, e.g., in 17.6% of subjects in, 33 but not in any model outcomes.An initially dominant Bacilli phase is associated in non-premature infants with a shorter gestational period, 33 but it is unclear exactly what factors are responsible.A similar initial dominance of Bacilli that often occurs in premature infants has been hypothesized to be related to selection pressures by the immune system, a different composition of initial colonizers, 39 or a defective mucin barrier. 40Secondly, the model predicted a very low Veillonella dispar abundance in all conditions.These predictions contradict in vivo data 21,41 in which V. dispar is relatively abundant.V. dispar likely has a lower abundance in the model due to an incorrectly reduced growth rate relative to the other species in the model on lactate, the primary energy source of V. dispar 42 (Figure 4B).We do not expect a large influence on the overall model predictions from this discrepancy, as C. acnes has a metabolism similar to that of V. dispar in the model and in vitro: both produce propionate, consume lactate, and cannot consume lactose. 38However, we cannot exclude that other species in the model, such as Veillonella spp., may be more important competitors in vivo than the competitors that the model predicts.
Despite the inevitable limitations of the model, we have shown here how the model can be used to produce testable predictions on the effects of prebiotics and competition on butyrate producing bacteria in the infant gut microbiota.Future versions of the model may be a useful help in follow-up studies on the effects of nutrition on bacterial population dynamics in the infant and adult gut microbiota.

Limitations of the study
Potential sources of the discrepancies between model predictions and experimental data include (1) errors in the metabolic predictions of the underlying FBA models; (2) computational errors, and (3) incomplete representation of the biology underlying infant digestion.A typical error occurring in FBA models is an incomplete prediction of metabolic shifts, which is in part due to the assumption of FBA models that the growth rate or energy production is optimized. 43For example, the FBA model does not correctly predict the metabolic shift from high-yield to lowyield metabolism as observed in vitro in Bifidobacterium growing on increasing concentrations of GOS and 2 0 -FL. 44,45FBA only predicts highyield metabolism.The model, therefore, likely underestimates total lactate production.The effects of this discrepancy on the results are difficult to predict, but as lactate is a cross-feeding substrate, the underestimation of lactate may cause the model to underestimate the abundance of cross-feeding species such as C. acnes or butyrate producing bacteria.The optimality assumption of FBA also ignores any other ''task'' that a bacterium has, besides growth.For example, sporulation, toxin production, or metabolic anticipation 46 may limit biomass production.The model does not represent such genetically regulated mechanisms.
Further errors in the model predictions can be due to simplifications in the FBA model.For example, we assume that the total flux through the reaction networks is capped by a limit that depends linearly on the local population size (Equation 5).This limit mimics the maximum volume in a cell that can be filled with enzymes.Here, each enzyme is assumed to have equal maximum flux, and the optimization problem then predicts the optimal relative flux distribution.In reality, due to differences in enzyme concentration and enzyme efficiency, these maximum fluxes can of course differ, which affects the predictions of FBA. 47,48If species-specific data on efficiency and genetic regulation of pathways become available, such weighting could be included in the model.The metabolic predictions from the FBA layer could be further improved in future versions of the model by improving the Gibbs free energy estimates (Equation 1), such as by using species-specific or dynamic estimates of intracellular pH.Intracellular pH is an important factor in calculating Gibbs free energy. 49Better Gibbs free energy estimates would then allow for reactions to be curated further based on whether they are thermodynamically feasible, and for thermodynamic favourability to be included in FBA.The inclusion of thermodynamic favourability in FBA has previously improved metabolic predictions for intracellular metabolism. 50,51Also, the FBA model currently includes an extracellular compartment in which long GOS chains are broken down into shorter GOS chains, but it is not possible for extracellular breakdown products to diffuse during this process.Such extracellular digestion may lead to additional competition effects, because competitors may ''steal'' digestion products without investing in the enzymes themselves. 52Such effects may become important if additional species are introduced in the model that digests prebiotics extracellularly, such as Bifidobacterium bifidum. 18omputational errors in the model (2) include the discretization of time, the discretization of space, and rounding errors in the FBA solver.Firstly, all processes in the model are assumed to be constant within each timestep, which means the model only roughly approximates the continuous temporal dynamics of processes such as metabolism and diffusion.Several techniques have been developed to improve this, such as incorporating additional rate of change constraints to the FBA. 53This ensures that fluxes can only change gradually between timesteps, thus bringing metabolism in the model closer to the continuous dynamics of actual bacterial metabolism.This may cause a mismatch between protein activity and the environment, where for example the proteins that take up lactose are still active even without lactose.Secondly, we discretize the three-dimensional continuous cylindrical space of the gut into a two-dimensional rectangular grid of lattice sites.We consider each lattice site to be of equal volume and to have equal flow through it.This simplification introduces many errors, as lattice sites must represent different shapes of three-dimensional space, and these shapes are not connected as they would be in three-dimensional space.It is difficult to estimate what impact these discretizations have on the model.Finally, the FBA solver uses floating point arithmetic to generate a deterministic but not exact solution to each FBA problem.These distortions are very small, typically on the order of 10 À 15 mmol per metabolite per FBA solution, so we do not expect a notable effect on the results.
Errors in the model predictions due to an incomplete representation of the biology underlying infant digestion (3) include missing organisms, missing ecological interactions, the simplifications we made to the metabolic input, and missing representation of host interactions.Firstly, the model does not include fungi or archaea in the infant gut.Both groups occur at a lower absolute abundance than the bacterial microbiota but may still influence it. 54Secondly, the model does not include interactions between bacteria other than cross-feeding and competition for resources.Missing interactions include acidification of the gut, 55 the production of bacteriocins, 56 and the effects of phage infections, 57 all of which have species-specific effects.Thirdly, the model does not include the input of fats, proteins, or minerals into the gut.This means that the model cannot represent stimulation of growth by digestion of fats or proteins, nor potential limits on growth due to, for example, the lack of iron 58 or essential amino acids. 59Finally, the model does not represent the interactions of the host with the microbiota, such as the continuous secretion by the gut wall of mucin 60 and oxygen, 61 and the uptake of short-chain fatty acids. 62Colonic mucins in particular could greatly influence the microbiota, as B. bifidum consumes colonic mucins extracellularly, which facilitates cross-feeding by butyrate producing bacteria in vitro. 63nulinivorans GEM.Roseburia spp.have been shown to be a prevalent butyrate producing bacterium in infants in other studies. 20Together, these form the list of species (Table 1).

Changes from AGORA
The model uses GEMs generated in the AGORA project. 25We have applied various changes and additions to these models (Table S1).
We have added digestion of GOS or 2 0 -FL to the B. longum ssp.infantis GEM as follows. 2 0 -FL digestion was implemented by adding reactions representing an ABC-transporter and an intracellular fucosidase that breaks 2 0 -FL down to lactose and fucose. 44GOS was represented through separate DP3, DP4, and DP5 fractions. 68The DP4 and DP5 fractions are broken down extracellularly to DP3 and DP4 fractions respectively, releasing one galactose molecule in the process. 69The DP3 fraction is taken up with an ABC transporter, and broken down internally to lactose and galactose. 69e have also further expanded earlier curation of the AGORA GEMs. 23We disabled anaerobic L-lactate uptake for the Bifidobacterium GEMs and for E. coli in line with available literature. 70,71To have the GEMs correspond with existing literature on lactose uptake we added a lactose symporter to Anaerobutyricum hallii, 63 both Bifidobacterium longum GEMs, 72 Roseburia inulinivorans, 73 Haemophilus parainfluenzae, 74 and Rothia mucilaginosa. 75We also added galactose metabolism to R. inulinivorans 76 and R. mucilaginosa. 75Further changes were made to prevent unrealistic growths and the destruction of atoms within reactions (Table S1).

Validity checks
After applying the changes in Table S1 we tested all GEMs individually for growth on a substrate of lactose and water.In line with literature, this did not lead to growth for Veillonella disparans, 42 Cutibacterium acnes, 38 Eggerthella sp.YY7918, 77 and Gemella morbillorum. 78All other species grew on this substrate.We also checked for any spurious growth by checking each GEM for growth with only water present.
During each simulation, the model checks the FBA solutions for thermodynamic plausibility.The model uses a database of Gibbs free energy values 49 for all metabolites except 2 0 -FL and GOS.Values for 2 0 -FL and GOS were generated from the values for lactose and fucose, and lactose and galactose, respectively.Separate values were generated for the separate fractions of GOS.All values assumed a pH of 7 and an ionic strength of 0.1 M.This is a simplification of actual pH values.The infant gut typically has a pH around 5-6.5. 79,80However, as most reactions in the model occur inside bacterial cells, we have used a pH of 7, which is a rough estimate of the internal pH of common gut bacteria at an external pH around 5-6.5. 81,82Ideally, pH would be calculated dynamically, and thermodynamic values would be generated dynamically based on the relevant pH values, but this has not been implemented in the current version of the model.
Energy loss l in joules per timestep per population unit is recorded as follows, where i are metabolites, F is the exchange flux rate in mol per timestep per population unit and E contains the Gibbs free energy in joules per mol for each metabolite, l = X i FðiÞ $ EðiÞ (Equation 1) We found that in the simulations of Figure 2A with the baseline level of lactose, combined with those with GOS and 2 0 -FL (n = 90) 99.98% of all FBA solutions had a lower or equal amount of Gibbs free energy in the output compared to the input.The remaining 0.02% of FBA solutions was responsible for 0.003% of total bacterial growth.

FBA with enzymatic constraint
Although other aspects of the model were changed, the FBA approach we used is identical to that used in the earlier model. 23The model uses a modified version of flux balance analysis with an enzymatic constraint to calculate the metabolic inputs and outputs of each population at each timestep. 29,43Each GEM is first converted to a stoichiometric matrix S. Reversible reactions are converted to two irreversible reactions, so that flux is always greater than or equal to 0. Reactions identified in the GEM as 'exchange', 'sink', or 'demand' in the GEM are also recorded as 'exchange' reactions.These exchange reactions are allowed to take up or deposit metabolites into the environment.Each timestep, all reactions are assumed to be in internal steady state: where f ! is a vector of the metabolic fluxes through each reaction in the network, in mol per time unit per population unit.
Each exchange reaction that takes up metabolites from the environment F in is constrained by an upper bound F ub which represents the availability of metabolites from the environment.It is determined as follows: where F ! in is a vector of fluxes between the environment and the bacterial population.F ! ub is a vector of upper bounds on these fluxes.F ! ub is set dynamically at each timestep t by the spatial environment at each lattice site x !: 4) where c ! is a vector of all metabolite concentrations in mol per lattice site, x ! is the location and Bð x ! ; tÞ is the size of the local bacterial population.The size of B can range from 5$10 7 to 2$10 10 bacterial cells.Finally the enzymatic constraint constrains the total flux through the network.It represents the maximum, total amount of flux that can be performed per cell in each population: X f !% a: (Equation 5) The enzymatic constraint a is in mol per time unit per population unit.As both f ! and a are per population unit, this limit scales linearly with population size, so each bacterial cell contributes equally to the metabolic flux possible in a lattice site.The enzymatic constraint is included as a constraint on each FBA solution.Given the constraints, FBA identifies the solution that optimizes the objective function, ATP production.
The solution consists of a set of input and output exchange fluxes F ! in ð x ! ; tÞ and F ! out ð x !;tÞ, and a growth rate gð x !;tÞ.The exchange fluxes are taken as the derivatives of a set of partial-differential equations to model the exchange of metabolites with the environment.The size of the population increases proportionally to the growth rate in the FBA solution.
To mimic quiescence at high densities, populations above the spreading threshold of 2$10 10 bacteria do not perform metabolism.In practice this rarely occurs because we maintain sufficient space for populations to spread into empty lattice sites.In the simulations of Figure 2A (n = 120) metabolism was not performed in, on average, 0.05% of all populations in a timestep.

Environmental metabolites
We model 735 different extracellular metabolites.This is the union of all metabolites that can be exchanged with the environment by at least one GEM in the model.In the simulations 39 metabolites are present in the medium in more than micromolar amounts at any point.We combine the L-lactate and D-lactate metabolites for Figure 1B and Videos S1 and S2.Nearly all lactate in the model is L-lactate.
To represent the mixing of metabolites by colonic contractions we apply a diffusion process to the metabolites at each timestep.Metabolic diffusion is applied in two equal steps to the model.In each step, 14.25% of each metabolite diffuses from each lattice site to each of the four nearest neighbors.This causes a net diffusion each timestep of 6:3$10 5 cm 2 /s.In the sensitivity analysis, eight of these diffusion steps are performed when the lattice is larger, or one step with a diffusion of 7.125% is performed when the lattice is smaller.This ensures the net diffusion of 6:3$10 5 cm 2 /s remains constant.When we use well-mixed conditions in the sensitivity analysis, all nutrients and metabolites are instead divided equally over all lattice sites.Metabolites are also added and removed by bacterial populations as a result of the FBA solutions, yielding where F ! out ð x ! ; tÞ is a vector of fluxes from the bacterial populations to the environment, in mol per time unit per population unit, D is the diffusion constant, L is the lattice side length, and NBð x !Þ are the four nearest neighbors.
All metabolites except oxygen are moved distally by one lattice site every timestep to represent advection.On the larger lattice the metabolites are instead moved by two lattice sites.On the smaller lattice, they are moved one lattice sites every two steps.The transit time, including diffusion, is approximately 11 h, corresponding with in vivo observations in newborn infants. 83,84Metabolites at the most distal column of the lattice, the end of the colon, are removed from the system at each timestep.
Every 60 timesteps (representing 3 h) metabolites representing inflow from the small intestine are inserted into the first six columns of lattice sites.Three hours represents a realistic feeding interval for neonates. 85Food intake contains 211 mmol of lactose by default, a concentration in line with human milk, 31 assuming 98% host uptake of carbohydrates before reaching the colon. 26In some simulations 211 mmol of additional lactose, GOS, or 2 0 -FL is added.Because there is very little uptake of prebiotics by the infant, 32 the oral intake of prebiotics would be much lower than that of lactose.GOS is inserted as separate fractions of DP3, DP4, or DP5 based on analysis of the composition of Vivinal-GOS. 6864% is DP3, 28% is DP4 and 8% is DP5.Water is provided in unlimited quantities.Oxygen is placed during initialization 86 at 0.1 mmol per lattice site.No other metabolites are available, other than those produced as a result of bacterial metabolism within the model.

Population dynamics
During initialization there is a probability of 0.3 for each lattice site to get a population of size 5$10 7 of a random species (Table 1).Taken together, this averages around 540 populations, leading to a total initial bacterial load of 2:7$10 10 , in line with in vivo estimates 87 when we assume a uniform bacterial density and a total colon volume of 90 mL.In each timestep each local population solves the FBA problem based on its own GEM, the enzymatic constraint a, its current population size Bð x ! ; tÞ and the local concentrations of metabolites c ! ð x !;tÞ, and applies the outcome to the environment (see above) and the growth rate gð x ! ; tÞ to its own population size, as follows: dBð x ! ; tÞ dt = Bð x ! ; tÞgð x ! ; tÞ: (Equation 7) Each step, each population of at least 1$10 10 bacteria (Table 2) will create a new population if an adjacent empty lattice site is available.Half of the old population size is transferred to the new population, so that the total size is preserved.To mimic colonization new populations are introduced at random into empty lattice sites during the simulation, representing both dormant bacteria from colonic crypts 88 and small bacterial populations formed from ingested bacteria, which may only become active after being moved far into the gut.Each empty lattice site has a probability of 0.00005 (Table 2) each step to acquire a new population of a randomly selected species.All species have an equal probability to be selected.We initialize these populations at the same population size B as the initial populations in the model (Table 2).Each population dies out at a probability of 0.0075 per timestep, creating a turnover within the range of estimated microbial turnover rates in the mouse microbiota. 89o mix the bacterial populations, the lattice sites swap population contents each timestep.We use an algorithm inspired by Kawasaki dynamics, 90 also used previously for bacterial mixing 23,24 In random order, the bacterial content of each site, i.e., the bacterial population represented by its size Bð x ! ; tÞ and the GEM, are swapped with a site randomly selected from the Moore neighborhood.This swap only occurs if both the origin and destination site have not already swapped in this timestep.With this mixing method the diffusion constant of the bacterial populations is 6:3$10 5 cm 2 =s, equal to that of the metabolites.Bacterial populations at the most distal column, i.e., at the exit of the colon, are removed from the system.To increase the bacterial diffusion rate in the sensitivity analysis this process was executed five times, marking all sites as unswapped after each execution.To decrease the bacterial diffusion rate the number of swaps was limited to a fifth of the usual number of swaps.In the sensitivity analysis simulations with a larger or smaller grid the bacterial diffusion rate was adapted in the same way, to ensure a constant effective diffusion in cm 2 =s.In the well-mixed sensitivity analysis simulations all populations were instead assigned random non-overlapping locations.

Data recording
We record the size, species, location, and important exchange fluxes F ! in ð x ! ; tÞ and F ! out ð x ! ; tÞ for each population at each timestep.To

Figure 1 .
Figure 1.Model predicts coexistence of Bifidobacterium and butyrate producing bacteria in absence of competition (A) Schematic of the model.Circles represent bacterial populations, color represents species.Flow through the tube is from left (proximal) to right (distal).Nutrients entered the system proximally.All metabolites leave the system distally.Lattice dimensions are schematic.(B) Screenshots of the model at a single time point, showing, from top to bottom, the bacterial layer, lactose, lactate, and acetate.Brightness indicates growth in the bacterial layer, and concentration in the metabolic layers.(Cand D) Abundance of (C) Bifidobacterium spp., (D) butyrate producing bacteria, at the end of 21 days for 30 sets of simulations with no prebiotics, no prebiotics and additional lactose, with GOS, or with 2'-FL at the end of 21 days n = 30 for each condition, each simulation is represented by one dot.See also TableS1.

Figure 2 .
Figure 2. Unlike GOS, 2 0 -FL leads to stimulation of butyrate producing bacteria through 1,2-PD in the full simulated microbiota For a Figure360 author presentation of this figure, see https://doi.org/10.1016/j.isci.2024.109085.(A) Relative abundance of bacterial species in the condition with no prebiotics, no prebiotics and additional lactose, with GOS, or with 2'-FL at the end of 21 days n = 30 for each condition, each simulation is weighed equally.The key to the species in each group is in Table 1.(B) Abundance of butyrate producing bacteria at the end of 21 days for the four conditions of A. n = 30 for each condition.Each simulation is represented by one dot.p < 0.001 for 2 0 -FL compared to no prebiotics and no prebiotics with additional lactose.p = 0.004 for 2 0 -FL compared to GOS. (C-E) Visualization of metabolic interactions in a sample simulation (C) without prebiotics (211 mmol lactose per 3 h) (D) with GOS (DP3,DP4, and DP5 displayed separately) (E) with 2 0 -FL.Line width is scaled with the flux per metabolite over the last 60 timesteps, multiplied by the carbon content of the molecule, with a minimum threshold of 100 mmol atomic carbon.Data from last 3 h, step 10,020 to 10,080.Circles indicate nutrients.(F)Abundance of butyrate producing bacteria with 2 0 FL at the end of 21 days.Uptake of lactose, lactate, or 1,2-PD by butyrate producing bacteria is disabled in the ''no lactose,'' ''no lactate,'' and ''no 1,2-PD'' conditions, respectively.p = 0.010,p < 0.001,p < 0.001 for each disabled uptake compared to the baseline, respectively.n = 30 for each condition.Each simulation is represented by one dot.NS: Not significant, *: p < 0.05, **:p < 0.01, ***:p < 0.001.See also FiguresS1 and S2and Videos S1, S2, S3, and S4.

Figure 3 . 2 0
Figure 3. 2 0 -FL makes butyrate producing bacteria resistant to competition by other infant gut bacteria (A) Abundance of butyrate producers with 2'-FL and without competitors (only Bifidobacterium and butyrate producers) at the end of 21 days.Uptake of lactose, lactate, or 1,2-PD is disabled for butyrate producers in the ''no lactose,'' ''no lactate,'' and ''no 1,2-PD'' conditions respectively.n = 30 for each condition.Each simulation is represented by one dot.(p < 0.001 for each disabled uptake compared to the baseline).(B-D) Abundance of butyrate producers at the end of 21 days (B) without prebiotics, either without competitors (only Bifidobacterium and butyrate producers), with addition of B. vulgatus, with addition of B. vulgatus unable to take up either lactose or lactate, or with addition of C. acnes.n = 30 for each condition.Each simulation is represented by one dot.p < 0.001 for abundance of butyrate producers with B. vulgatus compared to no competitors (C) with GOS, either without competitors (only Bifidobacterium and butyrate producers), with addition of C. acnes, with addition of C. acnes unable to take up lactate, or with addition of B. vulgatus.n = 30 for each condition.Each simulation is represented by one dot.p < 0.001 for abundance of butyrate producers with C. acnes compared to no competitors (D) with 2 0 -FL, either without competitors (only Bifidobacterium and butyrate producers), with addition of C. acnes, or with addition of B. vulgatus.n = 30 for each condition.Each simulation is represented by one dot.p = 0.001 for abundance of butyrate producers with C. acnes compared to no competitors.NS: Not significant, *: p < 0.05, **:p < 0.01, ***:p < 0.001.See also Figure S3.

Figure 4 .
Figure 4. Populations of butyrate producing bacteria only grow much faster than their competitors on a mixed substrate of 1,2-PD and lactate (A) Growth on unlimited lactose and water over a single timestep for butyrate producing bacteria (three rightmost bars, in green) compared to other lactosefermenting bacteria in the model.(B) Growth on unlimited lactate and water over a single timestep for butyrate producing bacteria (three rightmost bars, in green) compared to other lactatefermenting bacteria in the model.(C) Growth on unlimited 1,2-PD, acetate, and water over a single timestep for butyrate producing bacteria (two rightmost bars, in green) compared to another 1,2-PD-fermenting bacterial species in the model.(D) Growth on 1 mmol per mL of 1,2-PD and lactate, and unlimited acetate and water, over a single timestep for butyrate producing bacteria (three rightmost bars, in green) compared to other bacteria in the model for populations of 5$10 9 bacteria with access to one lattice site (0.05mL).

Table 1 .
Species and subspecies included in the model

Table 2 .
Parameters of the model